Cost - Effective Multiresolution Schemes for Shock Computations ∗

Harten’s Multiresolution framework has provided a fruitful environment for the development of adaptive codes for hyperbolic PDEs. The so-called cost-effective alternative [4, 8, 21] seeks to achieve savings in the computational cost of the underlying numerical technique, but not in the overall memory requirements of the code. Since the data structure of the basic algorithm does not need to be modified, it provides a set of tools that can be easily implemented into existing codes and that can be very useful in order to speed up the numerical simulations involved in the testing process that is associated to the development of new numerical schemes. In this paper we present two different applications of the cost-effective multilevel technique developed in [8] in non-standard situations. In both cases, its use has led to a drastic reduction in the computational time required for 2D numerical simulations, and hence to the ability to perform fine mesh simulation and, hence, to investigate new numerical techniques on personal computers. Introduction The numerical simulation of physical problems modeled by systems of conservation laws is rather delicate, due to the presence of discontinuities in the solution. Modern High-Resolution Shock Capturing (HRSC) schemes succeed in computing highly accurate numerical solutions, typically second, third or even higher order in smooth regions, while maintaining sharp, oscillation-free numerical profiles at discontinuities. HRSC schemes of the ENO-WENO family are well known nowadays in the hyperbolic community, due to their robust behavior in many different scenarios, even though no analytical proof of their convergence is available. There is a large body of scientific activity around these schemes, whose power lies in a sophisticated design of the numerical flux function in a conservative discrete formulation of the hyperbolic PDE or system. This is, in fact, the main drawback of these schemes, especially in multi-dimensional computations involving systems of PDEs, since the numerical flux computations tend to be fairly expensive. It is well known, however, that the costly numerical flux function of a HRSC scheme is only necessary in a neighborhood of singularities, hence a substantial amount of computational expense could be avoided by reducing the number of evaluations of costly numerical flux functions to those regions where singularities exist or are starting to develop. In the late eighties and early nineties (see e.g. [21]), Harten proposed various strategies to decide the regions in which a HRSC scheme must necessarily be used in order to maintain the high resolution properties in the numerical solution. Harten’s original multilevel strategies were developed for finite volume schemes for homogeneous conservation laws [4,21], where the numerical data are naturally treated ∗ Research supported in part by MTM2008-00974 1 Ecole Centrale Marseille and M2P2, UMR6181 2 Dept. Matematica Aplicada, Universitat de Valencia 3 Dept Matematica Aplicada I. Universidad de Sevilla c © EDP Sciences, SMAI 2009 Article published by EDP Sciences and available at http://www.edpsciences.org/proc or http://dx.doi.org/10.1051/proc/2009052 ESAIM: PROCEEDINGS 9 as approximations to the cell-averages of the true solution. Later on, Chiavassa and Donat [8] applied the cost-reduction strategy to Shu-Osher ENO-type finite-difference schemes, where the data are interpreted as point-values. This allows, in a rather natural way, to think of the basic ingredients of the smoothness analysis as interpolation errors, so that its relation to the smoothness of the underlying function is easy and well understood. The multilevel technique developed in [8], combined with a HRSC scheme developed in [13], has been applied in a number of situations with success. As observed in [8, 9, 32], the computational gain can reach a factor of 7, but normally its range lies between 3 and 5. In most situations, this is sufficient to obtain high resolution simulations on rather fine grids in 2D. The ability to obtain fine-mesh simulations in 2D is very relevant, since on many occasions diminishing the size of the computational mesh, and hence the associated numerical viscosity inherent to the scheme, can lead to numerical pathologies which are very hard to observe otherwise. In this work, we apply a straightforward extension of the multilevel strategy developed in [8] to nonhomogeneous systems of conservation laws. In particular, we shall investigate the properties of the extended multilevel technique on the shallow water system, and on a penalization model for compressible fluid flow with obstacles. We also mention here that there has been some recent work on the inclusion of source terms in fully adaptive multilevel strategies [24]. Fully adaptive schemes are not only cost-efficient, but also significant memory gains can be achieved through the full exploitation of the multilevel structure of the scheme, see [10] for a detailed account of this technique in the homogeneous case, and [33], [12] for later developments on these fully adaptive multilevel schemes. These techniques do require special data structures in order to obtain the expected memory gains, and its incorporation into an existing code is not straightforward. The cost effective alternative proposed in this work, can be incorporated almost as an external routine, at each time step, and remains the easiest multilevel scheme to adapt to an existing code. The rest of the paper is organized as follows: The basic cost-effective multilevel algorithm is briefly recalled in section 1. In section 2, we describe the shallow water system and the two numerical schemes we consider for its numerical approximation. In section 3 we present various numerical experiments that demonstrate that our multilevel strategy respects the well balancing properties of the underlying HRSC scheme, which is a fundamental property in shallow water simulations. In section 4 we recall a penalization technique for compressible flow with obstacles, developed in [5] and show several novel simulations that also illustrate the versatility of the numerical technique in obtaining fine-mesh simulations in non-standard situations. 1. The multilevel algorithm The general setting of the multilevel strategy we shall employ here has been described and analyzed in [8] for homogeneous conservation laws. In what follows we recall the main steps and explain our extension to the non-homogeneous case. Let us consider a 2D system of hyperbolic conservation laws ∂t~ U + ~ ∇ · ~ F(~ U) = B (1) where ~ U is the vector of conserved quantities, and B = B(~ U, ~ ∇~ U) represents the contribution of all terms not included in the divergence term. We shall only consider discretizations of this system on a Cartesian grid χ = {(xi = i4x, yj = j4y), i = 0, · · · , Nx j = 0, · · · , Ny} using the semi-discrete formulation ∂t~ Uij +Dij = Bij . (2) It is assumed that the numerical divergence is computed as Dij = ~ Fi+ 1 2 ,j − ~ Fi− 1 2 ,j 4x + ~ Gi,j+ 1 2 − ~ Gi,j− 1 2 4y , (3) 10 ESAIM: PROCEEDINGS where ~ F and ~ G are consistent numerical (normal) flux functions, which define the HRSC scheme used on a given application. The goal of the multilevel method is to reduce the CPU time associated to the underlying scheme by reducing the number of expensive flux evaluations involved in the computation of the numerical divergence. The computation of Bij depends on the problem at hand, and we shall present several situations in the following sections. By rewriting (2) as ∂t~ Uij + (Dij − Bij) = 0, (4) the basic mechanism of the cost-effective alternative for these problems is essentially that considered in [8], since the smoothness analysis of the solution vector ~ Uij reveals the behavior of the extended divergence Gij = Dij−Bij with respect to smoothness. 1.1. Smoothness analysis by Interpolatory Wavelet Transforms The use of an interpolatory wavelet transform in this context, as proposed in [8], is justified by the fact that in the finite-difference framework for ENO schemes proposed by Shu and Osher in [35], ~ Uij ≈ u(xi, yj). The Shu-Osher framework is used in all the numerical techniques we propose in this paper. In what follows, we shall restrict ourselves to 1D for simplicity in the description. The 2D extension is carried out by a standard tensor product construction. The different resolution levels are specified by a set of nested grids {χ, l = 0, 1, . . . , L}, where χ is considered the finest grid. The nested structure in the grid hierarchy implies xi ∈ χ ⇐⇒ x2li ∈ χ. If (v i )i are the point-values of a function v on χ, due to the embedding of the grids we have v i = v 0 2li, i = 0, . . . , Nl = Nx/2 , v i = v l−1 2i i ∈ Nl. (5) From the discrete data on χ a set of predicted values on χl−1 (the next finer level on the resolution ladder) can easily be computed by polynomial interpolation. Hence a set of predicted values, ṽl−1 i , can be computed by considering a centered interpolatory technique, in order to achieve maximal accuracy, ṽl−1 i = v l i/2 if xi ∈ χ l (6) ṽl−1 i = I [ xi; v ] = s ∑ k=1 βk ( v i+k−1 + v l i−k ) if xi ∈ χl−1\χl. The coefficients βk, k = 1, . . . , s can be easily determined (see e.g. [4]). The wavelet coefficients represent the difference between exact and predicted data, di = v l−1 i − ṽ l−1 i , xi ∈ χl−1 and, in this framework, they are interpolation errors which can be used directly as “regularity sensors” in order to localize non-smooth behavior. Notice that (6) implies that di = 0 for i ∈ χ, hence we can write vl−1 i = v l i/2, if xi ∈ χ l (7) vl−1 i = ṽ l−1 i + d l i if xi ∈ χl−1\χl. The sets {v} and {vl−1, d} are algebraically equivalent, and the same can be said about the sets v ↔ Mv = (v, d, . . . , d). Let us consider a Forward-Euler time-discretization of a 1D version of (4), i.e., U i = U n i −4tG i . (8) The application of any linear multiscale transformation M to this relation leads to MU i −MU n i = 4tMG i . (9) ESAIM: PROCEEDINGS 11 Hence, it is reasonable to assume that if a spatial location has associated small detail coefficients in both MU and MU, it would also have associated a small detail coefficient in MG. Hence, at this location, the value of G i might be correctly recovered simply by interpolating from lower resolution data. By using this heuristic argument, the information about the regularity of the data contained in MU is used to determine a flag vector (bi)l,i, whose value is either 0 or 1 for each spatial location according to the following rule: if |di| ≥ ε⇒ bi−k = 1 k,m = −2, . . . , 2 (10) if |di| ≥ 2ε and l > 0⇒ bl−1 2i−k = 1 k,m = −1, 0, 1 where r = 2s is the accuracy of the polynomial interpolatory technique used in the prediction step. The determination of the flag vector above follows Harten’s heuristics [21]. On one hand, it takes into account that large values of the detail coefficients correspond to non-smooth zones of the solution, like well resolved shocks or contact discontinuities. In addition, compression regions leading to shock formation, or steep profiles corresponding to discontinuities in the solution, exhibit a lack of regularity that can be estimated by looking at the decay rate of the detail coefficients across scales. Finally, a safety region of two cells around each flagged position is added. This produces a flagged region that should contain the singularities of both U and U (unknown at time n), since discontinuities of hyperbolic conservation laws move through the computational domain with a finite speed of propagation, and the CFL stability condition ensures that existing singularities do not travel farther than one cell per time step. 1.2. The multilevel evaluation of the discrete divergence The key point in the design of the multilevel algorithm in [8] consists in substituting the direct computation of the numerical divergence D on the finest grid by a multilevel strategy. The extended strategy we propose for the non-homogeneous case proceeds in the same manner, but instead of the values of the divergence operator, we apply the multilevel strategy proposed in [8] to the computation of the discrete G values. Notice that by using the smoothness analysis of U to compute the flag vector that determines the multilevel computation of G, we implicitly assume that the requirements of Harten’s heuristics remain unchanged: finite speed of propagation of singularities, controlled by the CFL condition for stability in the numerical scheme. The multilevel computation of G is, then, carried out as follows: Assuming that G is known on χ, the values of Gl−1 on χl−1 are computed as specified below • If xi ∈ χ, then Gl−1 i = G i/2. • If xi ∈ χl−1\χl, then Gl−1 i is computed using the boolean flag as follows: if bi = 1, compute Gl−1 i directly with the HRSC scheme if bi = 0, compute Dl−1 i = I [ xi;G ] = s ∑ k=1 βk ( G i+k−1 + G i−k ) . The process is repeated from l = L, . . . , 1 and, once it is completed, we obtain the values of G on the finest grid χ, which are needed by the ODE solver. 2. The shallow water equations The shallow water equations form a hyperbolic system of conservation laws that approximately describes various geophysical flows. When the effect of the bathymetry, or bottom surface elevation, is taken into account, the equations take the following form ∂t~ U + ∂x ~ F (~ U) + ∂y ~ G(~ U) = S(~x, ~ U) (11) 12 ESAIM: PROCEEDINGS where ~ U = (h, q1, q2) is the vector of unknowns, with h being the water depth and q1 = hu, q2 = hv the two components of the discharge ((u, v) is the velocity vector). The bottom topography is specified by the function z(x). The explicit expression of the 2D system is ∂t  h q1 q2 + ∂x  q1 q2 1 h + 12gh q1q2 h + ∂y  q2 q1q2 h q2 2 h + 1 2gh 2  =  0 −ghzx −ghzy  . (12) Non-trivial steady state solutions to (12) are physically relevant in many shallow water applications, and exist due to a balance of the flux divergence and the source term related to the variable bathymetry. Numerically preserving non-trivial steady states, or resolving small perturbations to them, when ∂x ~ F (~ U) + ∂y ~ G(~ U) ≈ S(~x, ~ U) yet both terms are relatively large, is a well known difficulty that has received considerable attention in the literature. In particular, the common fractional splitting approach to solve separately the convection and the source terms is known to produce spurious waves of a purely numerical nature. After many time steps, these disturbances can completely destroy the accuracy of the numerical solution. Much research has been devoted recently to finding numerical schemes that can accurately preserve non-trivial steady states at the discrete level. These schemes are called Well Balanced (WB) after the work of Greenberg and Le Roux [19, 20], or are said to satisfy the exact C-property, after the independent work of Bermudez and Vazquez-Cendon [3]. Well Balanced schemes, or schemes that preserve the C-property, incorporate an upwind treatment of the source term. This can be achieved by an appropriate de-centering of the source term contribution, as in e.g. [3], or by the use of specific Riemann solvers (e.g. [20, 25]) In this paper we shall consider two numerical schemes recently described in [6] and [27]. Both schemes incorporate an idea due to Gascon and Corberan [17] that becomes very useful in order to design well balanced schemes for balance laws. The balance law is first rewritten in a pseudo-conservative formulation, i.e. ut + f(u)x = s(x, u), ≡ ut + gx = 0 (13) with g(x, t) = f(u(x, t))− ∫ x x̄ s(y, u(y, t))dy. (14) Here x̄ is a reference point which can be taken as the initial point in the computational domain, for example. Both schemes follow the guidelines established by Shu and Osher in [35], the extension to 2D is done dimension by dimension, and the setting is specifically well adapted to the multilevel strategy described in the previous section. The specific details of both schemes are, however, quite different and, for the sake of completeness, we briefly review their essential design next. We refer the interested reader to [6, 27] for a full description of the schemes. 2.1. The 1J-2J scheme In [6] the authors extend the 2-Jacobian scheme developed in [13] for homogeneous systems of conservation laws (see also [14]) to the shallow water system. For high order simulations, the standard technique of separating the issues of spatial and temporal accuracy by a method of lines approach is followed, and the time discretization is carried out by a TVD-Runge Kutta scheme [35]. A brief algorithmic description of the numerical scheme is outlined next. ESAIM: PROCEEDINGS 13 For the scalar balance law (13), a method of lines discretization of the form ut + ĝi+ 1 2 − ĝi− 2 4x = 0 (15) is proposed, and the combined numerical fluxes ĝi± 2 evaluated in an ENO fashion, as specified in [35], using f ′(u) for determining the upwind direction. This procedure would involve integrations over domains of the form [x̄, xi] in the determination of the first-order terms in the expressions of ĝi± 2 , hence straightforward manipulations are carried out in order to arrive to an equivalent expression that only involves integration on consecutive cells. The equivalent form can be written as ut + g i+ 1 2 − g− i− 2 4x = 0. (16) It is observed in [6] that g i+ 2 − g− i+ 2 = ∫ xi+1 xi s(y, u(y, t))dy , (17) hence, the signed fluxes ĝ± i± 2 provide an upwind splitting of the source term contribution at the i+1/2 interface. The 1D shallow water system, Ut + F (U)x = S(U, x), (18) is first rewritten in the form Ut +G(U)x = 0, (19) with G(x, t) = F (U(x, t)) + B(x, t) and B(x, t) = (0, ∫ x x̄ ghzxds) . Then, the semi-discrete formulation of the numerical scheme for this 1D system can be expressed as Ut + G i+ 2 −G− i− 1 2 4x = 0, (20) where, again, the computation of G± i+ 2 only involves integral terms over consecutive cell centers. The construction of G± i+ 1 2 follows the basic design strategy of Marquina’s flux formula, developed in [13] for homogeneous conservation laws. At each cell-interface two states, U and U, and the corresponding spectral decompositions of the Jacobian matrix J(U) = ∂F/∂U , are computed at each side of the cell-interface by an appropriate nonlinear ENO interpolation process. Then, the numerical flux functions are obtained by applying the scalar algorithm to “sided” local characteristic variables and fluxes. Given U = U i+ 1 2 and U = U i+ 2 , the left and right states at the i+ 1 2 cell-interface, the flux functions G ± i+ 1 2 are computed as follows